Bound states of edge dislocations: The quantum dipole problem in two dimensions 
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We investigate bound state solutions of the 2D Schrodinger equation with a dipole potential 
originating from the elastic effects of a single edge dislocation. The knowledge of these states could 
be useful for understanding a wide variety of physical systems, including superfluid behavior along 
dislocations in solid 4 He. We present a review of the results obtained by previous workers together 
with an improved variational estimate of the ground state energy. We then numerically solve the 
eigenvalue problem and calculate the energy spectrum. In our dimensionless units, we find a ground 
state energy of -0.139, which is lower than any previous estimate. We also make successful contact 
with the behavior of the energy spectrum as derived from semiclassical considerations. 

PACS numbers: 67.80.B-, 02.60.-x 



I. INTRODUCTION 

There has been broad interest over the years in the 
physics of solids containing dislocations. In addition to 
affecting the mechanical properties of solids, the strain 
field associated with dislocations binds charge carriers in 
metals, or solute impurities in a generic solidP^ As such, 
the presence of dislocations has a significant effect on 
the transport, elastic and superconducting properties of 
the solid. In this context, it is important to know the 
spectrum of localized states due to a dislocation. 

In this article, we discuss the spectrum of bound states 
for an edge dislocation. Within linear elasticity theory 
the deformation potential due to an edge dislocation is 
proportional to the stress tensor or the divergence of the 
elastic displacement field. Considering a straight edge 
dislocation, oriented along the z-axis, within a continuum 
model, this potential is given by 



where p is the strength of the "dipole" potential, r is the 
distance from the dislocation axis and 9 is the azimuthal 
angle, both defined in the x-y plane. 1 The dipole mo- 
ment p depends on quantities such as the Fermi energy 
and the lattice and elastic constants of the solid. In an 
electrostatics context this potential can be realized as a 
dipole built by bringing two infinite line charges of op- 
posite sign close together. Here we address the quantum 
dipole problem by considering the solution of the corre- 
sponding two dimensional Schrodinger equation, 

h 2 „ 9 , cos# , _ . 

VV+P ib = Eib. (2) 

2m r 

For p > this potential is attractive for x < (thus, 
allowing for bound states) and repulsive for x > 0. It 
has parity in y; i.e., symmetry on reflection about the 
x-axis, which should be reflected in the cigenfunctions as 
well. The solution of the Schrodinger equation is com- 
plicated due to the non-central nature of the potential. 2 
The potential being non-separable further impairs the 
applicability of the WKB approximation. 



We are particularly motivated by the supersolid 
problem,^ and a possible interpretation of it which con- 
siders superfluidity to exist not in the bulk of solid 4 He 
but along a network of dislocations.^ We would like to 
solve the full nonlinear Ginzburg-Landau (GL) theory 
for such a system, for which we would first need to know 
the solution of the linearized equation. The lowest eigen- 
value of the linear equation is actually a measure of the 
local enhancement in T c produced by a dislocation. Fur- 
ther, the solution of the linear GL theory can be used 
to affect a separation of the transverse degrees of free- 
dom in the full nonlinear time-dependent GL equation, 
resulting in a one dimensional "amplitude equation" of 
superfluid density along the dislocation. The effective 
dimensionless coupling constant of this one dimensional 
theory is g = J dxdy\(f>o(x,y)\ , where (j>o(x,y) the nor- 
malized ground state eigenfunction of the linear GL equa- 
tion. Our numerical solution of the linear equation allows 
us to calculate this parameter which acts as an input to 
the weakly nonlinear analysis.^ 

The problem of finding the ground state energy of 
the quantum dipole problem has a long history, start- 
ing from the work of Landauer in 1954, 6 who used a 
variational approach. Subsequent authors used a va- 
riety of techniques for this estimate: semiclassicaP or 
purely variational^ methods, a combination of varia- 
tional and perturbative methods 10 or an expansion in 
terms of known basis functions] 11 * 12 ! but to our knowl- 
edge our work the first to solve t his u sing a direct nu- 
merical method. Some prior works^EH] have also studied 
the spectrum of the bound eigenstates. The ground state 
energies calculated in these works are shown in Table I, 
together with the numerical value obtained in this paper. 

In the next section, we provide the details of our vari- 
ational calculation for the ground state energy, including 
our choice of a suitable trial wave function. Next, we 
discuss the results and technical details of the several 
numerical methods we have used to calculate the eigen- 
value spectrum, and their relative merits and disadvan- 
tages. The methods involve diagonalization of the Hamil- 
tonian, carried out both in real space and in the basis of 
two dimensional hydrogenic wave functions (in contrast 
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TABLE I: Summary of ground state energy estimates of 
the edge dislocation potential. Energy is given in units of 
2mp 2 /h 2 . 





Ground state 
estimate 


Landauer (1954^ 


-0.102 


Emtage (1967^ 


-0.117 


Nabutovskii and Shapiro (1977^ 


-0.1014 


Slyusarev and Chishko (1984^ 


-0.1111 


Dubrovskii (1997^ 


-0.1196 


Farvacque and Francois (2001]P^ 


-0.1113 


Dorsey and Tonei^ 


-0.1199 


This work 


-0.139 



with previous calculations with different choices of basis 
expansions, e.g. Refs. [TT1 and IT2")) . Here we also compare 
the results obtained using different methods, which are 
found to agree in the essential features. In the final sec- 
tion, we provide a semiclassical argument to justify our 
results together with some discussion of the interesting 
properties of the classical problem. The semiclassical re- 
sult is found to extend to the lower energy eigenstates as 
well. 



II. VARIATIONAL CALCULATION 

Our initial approach to determine the ground state 
energy has been variational because this can be carried 
out analytically and provides a rough estimate which can 
then guide our more explicit numerical solution. Given a 
normalized wave function -0(r, 9), we minimized the en- 
ergy functional, 



d 2 x 



2m 



p- 



(3) 



This functional has its extrema at the solutions of the 
Schrodinger equation, Eq. Note that the length 

and energy scales which emerge from Eq. (J3l (or the 
Schrodinger equation) for this problem are a/2mp and 
2mp 2 /h 2 . In dimensionless variables, the normalized trial 
wave function used in our calculation is 



2AB (1 - r/BC) 




exp 



rcosflexp — — 

07T \ C 



(4) 



where A, B and C arc variational parameters. We choose 
the trial wave function so as to account for the anisotropy 
of the potential. Further, the asymptotic behavior of 
the potential is captured by the exponentially decaying 
factors. The minimum expectation value of the energy 
occurs when A = 0.803,5 = -0.774 and C = 2.14 with 



a ground state energy of —0.1199 which was found by 
Dorsey and TonerpSl This value is 2.5% lower than the 
previous lowest variational estimate (—0.1196) obtained 
by Dubrovskii PSl l n addition, by using this normalized 
trial wave function as the 4>o(x, y) we find the parameter 
g = j dxdy^^x.y)^ = 0.017. 



III. NUMERICAL METHODS 

A detailed numerical solution of the two-dimensional 
Schrodinger equation with the dipole potential, Eq. Q, 
is likely to provide more accurate ground state eigen- 
values in addition to determining the rest of the bound 
state eigenvalues and corresponding wave functions. We 
do this both by a real space diagonalization, where the 
Schrodinger equation is discretized on a square grid, and 
by expanding in the basis of the eigenfunctions of the 
two-dimensional Coulomb potential problem. Two spe- 
cial features of this dipole potential make it a numer- 
ically difficult problem: the singularity at the origin, 
and the long range behavior of the potential. It is ex- 
pected that the Coulomb wavefunctions would be better 
suited to capturing this long range behavior, and con- 
vergence would consequently be faster. Our results show 
that the Coulomb basis method is more accurate for the 
higher bound states (which are expected to extend more 
in space), as the real space methods are limited by size 
issues. However, the real space method works better for 
the ground state. 



A. Real Space Diagonalization Method 

For numerical purposes the Schrodinger equation is 
converted to a difference equation on a square grid of 
spacing h, with the Laplacian approximated by its five- 
point finite difference formJS] resulting in a block tridi- 
agonal matrix of size N 2 x N 2 , where the grid has di- 
mensions of TV x N. Each diagonal element corresponds 
to a grid point and has values of A/h 2 + V(x, y), whereas 
the nonzero offdiagonal elements all equal —1/h 2 . The 
matrix is thus very large but sparse. We use three 
different numerical methods to diagonalize this matrix: 
the biconjugate gradient method,^ the Jacobi-Davidson 
algorithm^ and Arnoldi-Lanczos algorithm,^ with the 
latter two being more suited to large sparse matri- 
ces whose extreme eigenvalues are required. We use 
freely available open source packages (JADAMILU 18 
and ARPACKP) written in FORTRAN for both. All 
three approaches are projective Krylov subspace meth- 
ods, which rely on repeated matrix- vector multiplications 
while searching for approximations to the required eigen- 
vector in a subspace of increasing dimensions. Refer- 
ence [20] provides a concise introduction to the Jacobi- 
Davidson method, together with comparisons to other 
similar methods. The implicitly restarted Arnoldi pack- 
age (ARPACK) is described in great detail in Ref. HU 
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Some general issues about the real space diagonalization 
as well as some specific features of the three methods 
used for it are discussed below. 

The accuracy of the real space diagonalization meth- 
ods is controlled by two main parameters: the grid spac- 
ing h and the total size of the grid, which is given by 
Nh. The finite difference approximation together with 
the rapid variation of the potential near the origin im- 
ply that the solution of the partial differential equation 
would be more accurate for a smaller grid spacing. We 
work with open boundary conditions, which means that 
a bound state wavefunction could be correctly captured 
only if the total size of the grid were to be greater than 
the natural decay length of the wavefunction. In other 
words, the eigenstate has to be given enough space to 
relax. This limits the number of bound states we can 
calculate accurately because a large grid size together 
with small grid spacings calls for a large number of grid 
points, thus quadratically increasing the size of the ma- 
trix to be diagonalized. Computational resources as well 
as the limitations of the algorithms themselves place an 
effective upper bound on the size of a diagonalizable ma- 
trix. We experimented to find that a 10 6 x 10 6 size sparse 
matrix was about the maximum that could be diagonal- 
ized with our computational resources. 

The origin of the square grid is symmetrically offset in 
both x and y directions to avoid the 1/r singularity. We 
tested first the accuracy of the real space techniques for 
the case of the two dimensional Coulomb potential, the 
spectrum of which is completely known. 22 We observe 
that for various lattice sizes the biconjugate method cap- 
tures at most the first four states whereas the Jacobi- 
Davidson method returns 20 eigenstates. The eigenval- 
ues obtained from both methods are accurate to within 
2% of the exact values.^ 

We have applied the biconjugate method to the edge 
dislocation potential for various lattice sizes, varying 
from 10x10 to 600x600. The number of eigenstates 
captured increases with the size of the lattice, as ex- 
pected. The ground state energy is observed to vary from 
-0.134 to -0.142. We also observe that for the number 
of grid points exceeding N = 2000 we encounter a nu- 
merical instability due to the accumulation of roundoff 
errors. For the largest real space grid size of 600x600 
(N = 1200, h = 0.5) we obtain seven eigenstates with a 
ground state energy of -0.1395. 

The ground state energy from the Jacobi-Davidson 
method, employed for the same lattice size gives -0.1395, 
which matches well with our expectations from the vari- 
ational calculation. We are able to obtain 20 bound state 
eigenvalues in this method using N = 1000, h = 0.5. It is 
checked that the low-lying eigenvalues are not very sen- 
sitive to values of h in this regime, so a relatively large 
value of 0.5 serves our purpose. 

The Arnoldi-Lanczos method yields very similar eigen- 
values. It takes more time and memory resources to con- 
verge but can calculate more eigenvalues, with greater 
accuracy for the higher excited states. It provides 30 
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FIG. 1: (Color online) Comparison of eigenvalues obtained 
from different methods. (The plot is on a log-log scale.) 



TABLE II: Comparison of first few energy eigenvalues ob- 
tained from different methods. Energy units: 2mp 2 /h 2 . n 
indicates quantum number of the state. 



n 


biconjugate 


Jacobi- 
Davidson 


Arnoldi- 
Lanczos 


Coulomb 
basis 


1 


-0.14 


-0.13954 


-0.13952 


-0.09697 


2 


-0.041 


-0.041480 


-0.041478 


-0.03281 


3 


-0.023 


-0.023314 


-0.023314 


-0.022067 


4 


-0.02 


-0.020086 


-0.020086 


-0.016744 


5 


-0.012 


-0.012592 


-0.012594 


-0.011944 



bound state eigenvalues for the same set of lattice param- 
eters as the above. Finally, after calculating the ground 
state wave function we find that the coupling constant 
g = 0.0194, slightly larger than the variational estimate 
of g = 0.017. 



B. Coulomb Basis Method 

We also calculate the spectrum numerically by using 
the linear variational method with the basis of the 2D hy- 
drogen atom wave functional There are two advantages 
of this method over the real space diagonalization meth- 
ods. First, the linear variational method is capable of 
capturing more excited states because the number of cal- 
culated bound states is not limited by the size of the real 
space grid but by the number of long-range basis func- 
tions. Second, the singularity at the origin of the edge 
dislocation potential docs not pose a problem anymore 
because elements of the Hamiltonian matrix become in- 
tegrable. 

Now we calculate the elements of the Hamiltonian 
matrix with a 2D edge dislocation potential. The 
Schrbdinger equation with the 2D Coulomb potential is 
analytically worked out in Ref. 1221 The normalized wave 
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functions of a 2D hydrogen atom are given by 

j— ( cos(Z0) for 1 < I < n, 

<iM) = ^Rn,i{r) x | ^ forZ=0, 

[sin(Z0) for -n < Z < -1, 

(5) 

where 



nA ' ~ (2\l\)\ V (2n - l)(n - |Z| - 1)! lPn J 

x exp (-^y) i F i(" n+ 1*1 + 1 ' 2 I / I + 1 " « r )> 

(6) 

with [3 n = 2/(2n — 1) and \F\ being the confluent hy- 
pergeometric function. The elements of the Hamiltonian 
with the 2D dipole potential are 



2 



4 J (7) 



1% 



cos 6 



7 / dr R niM (r)R n2 . h (r), 



(8) 

where = 5^ j a ±i/2 if both Zx and Z2 are less or greater 
than 0, or V — 1/V2 if Zi is and I2 positive or vice 
versa. The spectra are obtained for several total numbers 
of basis functions ./Vbasis- Due to the numerical precision 
in calculating elements of the Hamiltonian matrix ./Vbasis 
cannot be increased to more than 400. For -/Vb as is = 400 
we obtain about 149 bound states and the ground state 
energy of -0.0969. In order to improve the ground state 
energy, we introduce an additional decaying parameter 
in the basis functions, and optimize the energy levels for 
a certain value of this parameter. With the decaying 
parameter we obtain the best variational estimate for the 
ground state energy of -0.1257 for ./Vbasis = 400. 

We show the first twenty eigenvalues obtained from dif- 
ferent methods in FigJTIand the first five representative 
eigenvalues in Table |II| As mentioned earlier, the real 
space diagonalization methods provides a best estimate 
of the ground state energy whereas the Coulomb basis 
method is more suitable for higher excited states. The 
eigenvalues of both the Coulomb basis method and the 
real space diagonalization methods are found to match 
each other for excited states, and then they begin to de- 
viate again (see Fig. [2). This can be understood by the 
fact that the extent of wave functions of the 2D edge dis- 
location potential does not always increase as one goes to 
higher excited states - the wave functions of some excited 
states extend less than those of lower energy. There- 
fore, there are intermediate bound states that are to be 
missed in the real space calculation because the size of 
grid used in calculation is not large enough to capture 
them. For example, we find four more bound states with 
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FIG. 2: (Color online) Fit for the eigenvalue spectrum ob- 
tained from JADAMILU using f(x) — a(x — b) c . Fit values 
are —0.06, 0.61 and 0.96 for a, b and c respectively. 



the Coulomb basis calculation between the 18 th and 19 th 
excited states as calculated from the real space diagonal- 
ization method. This feature also explains the abrupt 
increase of the eigenvalue of the 19 th state calculated by 
using the Arnoldi-Lanczos method (ARPACK routine) in 

Fig-H 



IV. SEMICLASSICAL ANALYSIS 

It is usually insightful to consider the semiclassical so- 
lution of a quantum mechanics problem, since the higher 
energy eigenstates tend to approach classical behavior. A 
semiclassical estimate of the energy spectrum has been 
provided in Ref. [71 Here the total number of eigenstates 
up to a value of energy E is proportional to the volume 
occupied by the system in the classical phase space. This 
is expressed by Weyl's theorenPH ; 



n(E) 



A 2m\E\ 

4^ n 2 



o 



2mp 2 



\E\ 



(9) 



where A is the classically accessible area in real space 
and \E\ the absolute value of energy of the state. The 
higher order corrections can be shown to be less impor- 
tant for higher excited states, which is where the semi- 
classical picture applies. To find A, we need the classical 
turning points for this potential determined by setting 
E = V(r,6). Then the accessible area is the interior 
of a circle given by (x — ^e) 2 + y 2 — (^e) 2 , with area 
A = ir(p/2E) 2 . Therefore, we obtain (writing the nondi- 
mensionalized energy in our system of units as e) : 



n(e) = 



1 

16? 



(10) 



where n is the quantum number of the eigenstate, and e 
the corresponding energy. Note that the density of states 
dn/de scales as 1/e 2 . 
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wavefunctions obtained from our numerical calculations 
have been included in Fig. [3j The parity of the potential 
shows up in the wavefunctions being either symmetric or 
antisymmetric about the £-axis, although states of such 
"odd" and "even" parity do not always alternate. Also, 
as mentioned earlier, the spatial extent of the wavefunc- 
tions does not scale monotonically with quantum num- 
ber. We do not have any satisfactory explanation yet for 
these irregular features. 



(a)Eigenfunction for the ground state 
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(b)Eigenfunction for the 1st excited 
state 



FIG. 3: (Color online) Eigenfunctions of the first two bound 
states 



To check this result we fit the numerical spectrum with 
the following functional form: 



e(n) — a(n — b) c 



(11) 



with the fitting parameters having values a = —0.06, b — 
0.5, c — —0.98, each correct to within 5%. (Since we are 
dealing with bound states here, all the energy eigenval- 
ues are negative, and the higher excited states have lower 
absolute eigenvalues.) We show the fit to the spectrum 
obtained from JADAMILU routine in Fig. [2j The semi- 
classically derived dependence is found to closely match 
with the fit for numerically calculated energy cigcnstatcs, 
except for the b — 0.5 factor. In the limit of large n val- 
ues i.e, higher excited states, the fit relation tends to the 
scmiclassical result as expected. 

The classical trajectories for this potential bear the sig- 
nature of chaotic dynamics showing space-filling nature 
and strong dependence on initial conditions. However, 
for reasons not yet clear to us, they are not ergodic, filling 
up only a wedge-shaped region in real space instead of the 
full classically allowed circle. The quantum mechanical 
probability density as calculated from the eigenfunctions 
also exhibits such wedge-shaped regions. Some sample 



V. SUMMARY 

In conclusion, we have investigated the longstanding 
quantum problem of a two-dimensional dipole potential. 
The wave functions and the spectrum are calculated by 
solving the Schrodinger equation with the 2D dipole po- 
tential numerically, and also, in the case of the ground 
state, variationally. We find that the results obtained 
from the different methods are consistent and compare 
favorably with previous estimates in the literature. We 
also discover a simple pattern in the spectrum, (n cx e _1 ), 
which can be justified from scmiclassical considerations. 
Certain features of the spectrum and wave functions are 
yet to be explained and might provide scope for future in- 
vestigation. For example, the statistics of the level spac- 
ings could possibly be a signature of quantum chaos. We 
hope to extend our work to studying dislocation-induced 
superfluidity as a model of 4 He supersolidP In such a 
model, the linearized GL equation is isomorphic to the 
Schrodinger equation, Eq. ([2]), and the ground state en- 
ergy and its wave function determined in this work pro- 
vides an input to obtain the coupling constant g and the 
superfluid transition temperature modified by the pres- 
ence of dislocations. 
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